Importando as bibliotecas que serão usadas e definindo pi = π
import math
import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
import scipy.io as spio
import IPython.display as ipd
pi = np.pi
Ignorando os warnings causados pelas funções kaiser() e wavfile.read()
import warnings
warnings.filterwarnings("ignore", category=SyntaxWarning)
warnings.filterwarnings("ignore", category=spio.wavfile.WavFileWarning)
Copiando a função kaiser()
def kaiser(wp,wr):
wc = (wp + wr)/2
d = 0.01
Ap = 20*math.log10((1+d)/(1-d))
Ar = -20*math.log10(d)
if Ar < 21:
b = 0
D = .9222
elif Ar < 50:
b = 0.5842*(Ar-21)**0.4+0.07886*(Ar-21)
D = (Ar - 7.95)/14.36
else:
b = .1102*(Ar-8.7)
D = (Ar - 7.95)/14.36
k = math.ceil(math.pi*D/(wr-wp)-.5)
M = 2*k+1
n = np.arange(-k,k+1,1)
w = np.i0(b*np.sqrt(1-(4/M**2)*n**2))
w = np.divide(w,np.i0(b))
h = wc/math.pi*np.sinc(wc*n/math.pi)*w
return h
Importando a função espectro() (PS: Essa função será chamada com ';' no final para suprimir o output indesejado)
""" Rotina que exibe o espectro de magnitude (X(ejw)) de um sinal discreto """
def espectro(y):
#modulo da transf. de Fourier
Y = np.abs(np.fft.fft(y))
#frequencias avaliadas
w = np.linspace(0,2*math.pi,Y.size)
#exibe o grafico do espectro
plt.figure()
plt.plot(w,Y)
plt.xlabel('$\Omega$ [rad]', fontsize=10)
plt.ylabel('|$Y(e^{j\Omega})$|', fontsize=10)
plt.grid(True)
plt.xlim((0,2*math.pi))
plt.show()
return Y,w
1) Amostragem e Aliasing¶
Carregando o arquivo de áudio e tornando-o Mono
Fs, y = spio.wavfile.read("creed_overcome.wav")
y = (y[:,0] + y[:,1]) / 2
a) Análise do conteúdo espectral¶
Utilizando a função espectro para plotar o áudio em um gráfico
espectro(y);
Discussão:¶
Nota-se que há uma maior concentração de energia nos valores de Ω próximos de múltiplos de 2π, indicando que a maior parte da energia está nas frequências graves.
De mesmo modo, nota-se uma menor concentração de energia nos valores de Ω próximo de π, indicando que os componentes de alta frequência, ou seja, frequências mais agudas, têm menores intensidades.
Sobre a relação da frequência digital Ω com a frequência analógica $f$, sabe-se que: $$f = \frac{\Omega}{2\pi} \cdot F_s$$ Com $F_s$ sendo a taxa de amostragem.
b) Decimação do sinal com M = 6¶
Implementando uma função que faz a decimação
def decimacao(y, m):
return y[::m]
Fazendo a decimação com M = 6
y_dec = decimacao(y, 6)
Plotando o espectro do sinal decimado y_dec
espectro(y_dec);
Discussão:¶
Nota-se que a energia apresenta-se de maneira mais bem distribuída ao longo do gráfico. Os pontos próximos de π não estão mais tão abaixo que os pontos próximos de 0 e 2π. Isso é resultado do aliasing e denota a perda de alguns aspectos importantes da música original, como alguns agudos.
Nota-se também que a amplitude de todos os pontos está bem menor pois, ao decimar com fator M = 6, há uma divisão dessa amplitude pelo mesmo fator, uma vez que a energia total do sinal também é reduzida nessa proporção. O que é evidente uma vez que os valores mais altos do espectro do sinal original estavam em torno de 24000 e agora estão em torno de 4000. Vale ressaltar que esse diminuição não representa uma perda de energia real, apenas uma escala normalizada adequada ao novo domínio decimado.
Fazendo o display dos áudios y e y_dec, respectivamente
ipd.display(ipd.Audio(y, rate=Fs))
Fs_dec = Fs / 6
ipd.display(ipd.Audio(y_dec, rate=Fs_dec))
Discussão:¶
É visível como o áudio decimado perde qualidade em diversos aspectos. Conforme previsto, há uma tendência de "homogenização" do áudio, uma vez que ele perde aspectos mais agudos. Um exemplo notável são os pratos da bateria que, no original, são marcantes e constantes, mas no decimado se mostram ausentes. Após a decimação, o som se torna mais abafado, distante, áspero e metálico.
Importante ressaltar também a presença de uma espécie de chiado ou ruído quando o vocalista canta palavras com "s". Isso é bem notável próximo da marca de tempo 0:21, durante o trecho "now it's my turn to speak".
c) Minimização do Aliasing¶
Criando FPB para Ωp = 0,45 [rad] e Ωr = 2 [rad]
h1 = kaiser(0.45, 2)
Criando FPB para Ωp = 0,45 [rad] e Ωr = 0,50 [rad]
h2 = kaiser(0.45, 0.50)
Criando FPB para Ωp = 1,50 [rad] e Ωr = 2 [rad]
h3 = kaiser(1.50, 2)
Plotando os filtros em um gráfico
espectro(h1)
espectro(h2)
espectro(h3);
Discussão:¶
Nota-se, claramente, que o FPB gerado pelos valores Ωp = 0,45 [rad] e Ωr = 0,50 [rad] é o mais próximo de um FPB ideal, uma vez que este apresenta pouca oscilação antes e depois da banda de rejeição, bem como apresenta uma queda suficientemente brusca de 0 para 1.
d) Filtrando o sinal¶
Convoluindo o sinal y com o FPB h2, gerado com os valores Ωp = 0,45 [rad] e Ωr = 0,50 [rad]
saida = sp.signal.fftconvolve(y, h2, mode="full")
Plotando o espectro da saída após a convolução com o FPB
espectro(saida);
Discussão:¶
Nota-se que a concentração da energia nos pontos próximos de 0 e 2π se tornou ainda maior, uma vez que, após passar pelo FPB, conforme o esperado, as frequências altas próximas de π foram barradas, atenunando ainda mais as frequências baixas nos cantos do gráfico.
Toda a região do meio do gráfico apresenta-se basicamente nula. Isso causa a mesma perda de agudos que a decimação feita no item c) causou, mas prepara o terreno para que o aliasing não ocorra, uma vez que já espurga as frequências altas agora.
Fazendo o display da saida
ipd.display(ipd.Audio(saida, rate=Fs))
Discussão:¶
Nota-se que, conforme esperado, os sons agudos, como os pratos da bateria, não são possíveis de se escutar. A qualidade do áudio caiu, tornando-o mais distante e vazio, e alguns sons agudos se perderam, mas a voz do vocalista e os instrumentos mais graves não apresentaram tanta distorção.
e) Decimando a saída após passar pelo FPB¶
Decimando o sinal saida para M = 6
saida_dec = decimacao(saida, 6)
Plotando o espectro da saida_dec e do y_dec
espectro(saida_dec)
espectro(y_dec);
Discussão:¶
Nota-se que o espectro do sinal saida_dec, nos valores perto de π, se aproximam bem mais de 0, indicando uma supressão maior das altas frequências, mesmo após a decimação, ao contrário do espectro do sinal y_dec, que por sua vez apresenta uma estrutura bem mais homogênea perto de π.
Essa diferença nos valores perto de π gera diferenças, também, nos valores próximos de 0 e 2π. Na saida_dec, os picos se aproximam de 5000, enquanto na y_dec os picos não chegam nem a 4000.
Essas diferenças, em tese, resultarão em uma redução do aliasing, bem como uma maior preservação da estrutura do sinal original.
Fazendo o display da y_dec e da saida_dec
ipd.display(ipd.Audio(y_dec, rate=Fs_dec))
ipd.display(ipd.Audio(saida_dec, rate=Fs_dec))
Discussão:¶
Escutando os dois áudios lado a lado, é evidente que, conforme esperado, os agudos foram perdidos em ambos (pratos da bateria), mas em saida_dec, ao contrário de y_dec, não há a distorção causada pelo aliasing nos momentos em que o vocalista canta palavras com "s". Isso é visível, novamente, na marcação de tempo 0:21, pois no trecho "now it's my time to speak" o chiado ou ruído não acontece mais.
Mesmo que saida_dec também perca os agudos e, consequentemente, faça com que o som fique mais distante e com menos instrumentos, a ausência do aliasing torna o áudio mais bem preservado e agradável de se escutar que o y_dec.
2) DFT e Identificação de Frequências¶
f) Análise de um EEG¶
Abrindo o arquivo .csv e definindo a taxa de amostragem Fs
eeg = np.loadtxt("EEG.txt", dtype=float)
Fs = 250
Calculando a transformada de Fourier de eeg
n = len(eeg)
eeg_fft = np.fft.fft(eeg)
eeg_fft_mag = np.abs(eeg_fft) / n
Calculando o vetor de frequências
freqs = np.fft.fftfreq(n, 1/Fs)
Considerando apenas o espectro positivo de eeg_fft_mag e freqs
half = n // 2
eeg_fft_mag_half = 2 * eeg_fft_mag[:half]
freqs_half = freqs[:half]
Plotando o EEG num gráfico
plt.plot(freqs_half, eeg_fft_mag_half)
plt.xlabel("Frequência Analógica (Hz)")
plt.ylabel("Magnitude")
plt.title("Espectro de Magnitude do EEG")
plt.grid(True, alpha=0.3)
plt.show()
Indentificando a frequência de maior pico
k_max = np.argmax(eeg_fft_mag_half[1:]) + 1
Fe = freqs_half[k_max]
print(f"Frequência de maior pico = {Fe:.2f} Hz")
Frequência de maior pico = 11.67 Hz